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Quasiparticles in the vortex lattice of strongly type-II superconductors are investigated by means of 
a singular gauge transformation applied to the tight binding lattice Bogoliubov-de Gennes Hamil- 
tonian. We present a detailed derivation of the gauge invariant effective low energy Ifamiltonian for 
the quasiparticle- vortex system and show how the physics of the "Doppler shift" and "Berry phase" 
can be incorporated at the Hamiltonian level by working in the singular gauge. In particular, we 
show that the "Berry phase" effect manifests itself in the effective Hamiltonian through a half-flux 
Aharonov-Bohm scattering of quasiparticles off vortices and stress the important role that this 
effect plays in the quasiparticle dynamics. Full numerical solutions in the regime of intermediate 
fields Hc\ <C -B <C Hc2 are presented for model superconductors with s-, p- and d-wave symme- 
tries and with square and triangular vortex lattices. For s- and p-wave cases we obtain low energy 
bound states in the core, in agreement with the existing results. For d-wave case only extended 
quasiparticle states exist. We investigate in detail the nature of these extended states and provide 
comparison to the previous results within linearized "Dirac fermion" model. We also investigate 
internodal interference effects when vortex and ionic lattices have high degree of commensurability 
and discuss various specific choices for the singular gauge transformation. 



I. INTRODUCTION 

In conventional (s-wave) superconductors the single 
particle fermionic excitations (quasiparticles) are fully 
gapped everywhere on the Fermi surface and the quasi- 
particle density of states vanishes below a specific en- 
ergy. This has profound consequences for the tradi- 
tional phenomenology of superconductors. The gap in 
the fermionic spectrum leads to the well known activated 
BCS form of the quasiparticle contribution to various 
thermodynamic and transport properties. Furthermore, 
even as one moves beyond the mean-field BCS theory, the 
absence of low-energy quasiparticles in the superconduct- 
ing state allows one to rewrite the problem of supercon- 
ducting fluctuations as a "bosonic" theory, with the role 
of bosons played by fluctuating Cooper pairs, after inte- 
grating out "fermionic" degrees of freedom, i.e. the quasi- 
particles. In high temperature superconductors (HTS), 
however, everything is different: the cuprates appear to 
be accurately described by the dx2-y2-weLve order param- 
eter 0, consequently allowing quasiparticle excitations 
at arbitrary low energy near the nodal points. These 
low energy fermionic excitations appear to govern much 
of the thermodynamics and transport in the HTS mate- 
rials. We are thus handed a new intellectual challenge 
1^: we must devise methods that can incorporate the 
low energy fermionic excitations into the phenomenology 
of superconductors, both within the mean-fleld BCS-like 
theory and beyond. 

This challenge is not trivial and has many diverse com- 
ponents: low energy quasiparticles are scattered by im- 
purities in novel and unusual ways, depending on the 
low energy density of states they interact with ex- 
ternal perturbations in ways not encountered in conven- 
tional superconductors and these interactions give rise 
to new phenomena the low energy quasiparticles 



are expected to qualitatively affect the quantum critical 
behavior of HTS. Among many aspects of this new quasi- 
particle phenomenology a particularly prominent role is 
played by the low lying quasiparticle excitations in the 
mixed (or vortex) state. All HTS are extreme type-II 
systems and have a huge mixed phase extending from 
the lower critical field Hd which is in the range of 10- 
100 Gauss to the upper critical field fl"c2 which can be as 
large as 100-200 T. We suspect that in this large region 
the interactions between quasiparticles and vortices play 
the essential role in defining the nature of thermodynamic 
and transport properties. 

Such thermodynamic and transport properties are ex- 
pected to be rather different for distinct classes of uncon- 
ventional superconductors. This difference stems from 
a complex motion of the quasiparticles under the com- 
bined effects of both the magnetic field B and the local 
drift produced by chiral supercurrents of the vortex state. 
For example, in HTS the d^2_y2--wave nature of the gap 
function results in its vanishing along nodal directions. 
Along these nodal directions the pair-breaking induced 
by supercurrents has a particularly strong effect. On the 
other hand, in unconventional superconductors with the 
Px±iy pairing, Sr2Ru04 being a possible candidate 
the spectrum is fully gapped but the order parameter is 
chiral even in the absence of external magnetic field. This 
leads to two different types of vortices for two different 
field orientations |@,||. 

Still, in all these different situations, the quantum dy- 
namics of quasiparticles in the vortex state contains two 
essential common ingredients. First, there is a purely 
classical effect of a Doppler shift [Q,D: a quasiparti- 
cle energy is shifted by a locally drifting superfluid, 
i?(k) £^(k) — tiVs{r) ■ k, where Vs(r) is the local su- 
perfluid velocity. Vs(r) contains information about vor- 
tex configurations allowing us to connect quasiparticle 
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spectral properties to various cooperative phenomena in 
the system of vortices [p|-|lT|. The Doppler shift effect 
is not pecuhar to the vortex state. It also occurs in the 
Meissner phase and is generally present whenever a 
quasiparticle experiences a locally uniform drift in the 
superfluid velocity. Second, there is also a purely quan- 
tum effect which is intimately tied to the vortex state: as 
a quasiparticle circles around a vortex while maintaining 
its quantum coherence, the accumulated phase through a 
Doppler shift is ±7r. This implies that there must be an 
additional compensating ±7r contribution to the phase 
on top of the one due to the Doppler shift [|l2|. The 
required ±7r contribution is supplied by a "statistical in- 
teraction" or a "Berry phase" effect and can be built in at 
the Hamiltonian level as a half-flux Aharonov-Bohm scat- 
tering of quasiparticles by vortices | [T^ ]. This interplay 
between the classical (Doppler shift) and purely quan- 
tum effect ( "Berry phase" ) is what makes the problem of 
quasiparticle- vortex interaction particularly fascinating. 

Let us briefly review what is already known about the 
subject. The initial theoretical investigations of gapped 
and gapless superconductors in the vortex state were di- 
rected along rather separate lines. The low energy quasi- 
particle spectrum of an s-wave superconductor in the 
mixed state was originally studied by Caroli, de Gennes 
and Matricon (CdGM) ||l^ within the framework of the 
Bogoliubov-de Gennes equations Their solution 

yields well known bound states in the vortex cores. These 
states are localized in the core and have an exponential 
envelope the scale of which is set by the BCS coherence 
length. The low energy end of the spectrum is given by 

~ fi{AQ/Ep), where fi = 1/2,3/2,..., Aq is the over- 
all BCS gap and Ep is the Fermi energy. This solution 
can be relatively straightforwardly generalized to a fully 
gapped, chiral p-wave superconductor. In this case the 
low energy quasiparticle spectrum also displays bound 
vortex core states, whose energy quantization is, however, 
modified relative to its s-wave counterpart, precisely be- 
cause of the chiral character of a px±iy-'wai,ve supercon- 
ductor and the ensuing shift in the angular momentum. 
For example, the low energy spectrum of quasiparticles 
in the singly quantized vortex of the px±iy-'wave super- 
conductor, possesses a state at exactly zero energy 0,1). 

By comparison, the spectrum of a gapless d-wave su- 
perconductor in the mixed phase has become the subject 
of an active debate only relatively recently, fueled by the 
interest in HTS. Naturally, the first question that arises 
is what is the analog of the CdGM solution for a single 
vortex. It is important to realize here that the situation 
in a dj.2_y2 superconductor is qualitatively different from 
the classic s-wave case : when the pairing state has a 
finite angular momentum and is not a global eigenstate 
of the angular momentum (a dx^-y-i superconductor 
is an equal admixture of = ±1 states), the problem 
of fermionic excitations in the core cannot be reduced 
to a collection of decoupled ID dimensional eigenvalue 
equations for each angular momentum channel, the key 
feature of the CdGM solution. Instead, all channels re- 



main coupled and one must solve a full 2D problem. The 
fully self-consistent numerical solution of the BdG equa- 
tions 16 1 reveals the most important physical conse- 
quence of this qualitatively new situation: the vortex core 
quasiparticle states in a pure dx2_y2 superconductor are 
delocalized with wave functions extended along the nodal 
directions. The low lying states have a continuous spec- 
trum and, in a broad range of parameters, do not seem 
to exhibit strong resonant behavior. Obviously, this is in 
sharp contrast with a discrete spectrum and true bound 
quasiparticle states of the CdGM s-wave solution. We 
expect the above qualitative results to hold for all un- 
conventional superconductors and within confines of the 
simple BdG equations, as long as there are nodes in the 
gap. 

A particularly important issue in this context is the na- 
ture of the quasiparticle excitations at very low fields, in 
the presence of a vortex lattice. This is a novel challenge 
since the spectrum starts as gapless at zero field and at 
issue is the interaction of these low lying quasiparticles 
with the vortex lattice. This problem has been addressed 
via numerical solution of the tight binding model |l7 , a 
numerical diagonalization of the continuum model IS] 
and a semiclassical analysis Q. Gorkov, Schrieffer 1£] 
and, in a somewhat different context, Anderson pre- 
dicted that the quasiparticle spectrum is described by a 
Dirac-like Landau quantization of energy levels 



En = ±hujH\/n, 



0,1,. 



(1) 



where ujh — ^/^LO^^KQ/h, Uc = eB/mc is the cyclotron 
frequency and Aq is the maximum superconducting gap. 
The Dirac-like spectrum of Landau levels arises from 
the linear dispersion of nodal quasiparticles at zero field. 
This argumentation neglects the effect of spatially vary- 
ing supercurrents in the vortex array which were shown 
to strongly mix individual Landau levels [pi[ . 

Recently, Franz and Tesanovic (FT) |12f| pointed out 
that the low energy quasiparticle states of a (i^2_j,2- 
wave superconductor in a vortex state are most natu- 
rally described by strongly dispersive Bloch waves. This 
conclusion was based on the particular choice of a sin- 
gular gauge transformation, which allows for the treat- 
ment of the uniform external magnetic field and the ef- 
fects produced by chiral supercurrents on equal footing. 
The starting point was the Bogoliubov-de Gennes (BdG) 
equation linearized around a Dirac node. By employ- 
ing the singular gauge transformation FT mapped the 
original problem onto that of a Dirac Hamiltonian in pe- 
riodic vector and scalar potentials, comprised of an ar- 
ray of an effective magnetic Aharonov-Bohm half-fluxes, 
and with a vanishing overall magnetic flux per unit cell. 
The FT gauge transformation allows use of standard 
band structure and other zero-field techniques to study 
the quasiparticle dynamics in the presence of vortex ar- 
rays, ordered or disordered. Its utility was illustrated in 
Ref. JT^ through computation of the quasiparticle spec- 
tra of a square vortex lattice. A remarkable feature of 
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these spectra is the persistence of the massless Dirac 
node at finite fields and the appearance of the "hnes of 
nodes" in the gap at large values of the anisotropy ratio 
ao ~ vp/v/^, starting at au — 15. Furthermore, the FT 
transformation directly reveals that a quasiparticle mov- 
ing coherently through a vortex array experiences not 
only a Doppler shift caused by circulating supercurrents 
but also an additional^ "Berry phase" effect: the latter 
is a purely quantum mechanical phenomenon and is ab- 
sent from a typical semiclassical approach. Interestingly, 
the cyclotron motion in Dirac cones is entirely caused 
by such "Berry phase" effect, which takes the form of a 
half-flux Aharonov-Bohm scattering of quasiparticles by 
vortices, and does nof explicitly involve the external mag- 
netic field. It is for this reason that the Dirac-like Landau 
level quantization is absent from the exact quasiparticle 
spectrum. 

Further progress was achieved by Marinelli, Halperin 
and Simon |£2| who presented a detailed perturbative 
analysis of the linearized Hamiltonian of Ref . . They 
showed that the presence of the particle-hole symmetry 
is of key importance in determining the nature of the 
spectrum of low energy excitations. If the vortices are ar- 
ranged in a Bravais lattice, they showed that, to all orders 
in perturbation theory, the Dirac node is preserved at fi- 
nite fields, i.e the quasiparticle spectrum remains gapless 
at the r point. This result masks intense mixing of in- 
dividual basis vectors (in the case of Rcf. these are 
Dirac plane waves) , including strong mixing of states far 
removed in energy. The continuing presence of the mass- 
less Dirac node at the T point after the application of 
the external field is thus not due to the lack of scatter- 
ing which is actually remarkably strong. Rather, it is 
dictated by symmetry: Marinelli et al. demonstrated 
that the crucial agent responsible for the presence of the 
Dirac node is the particle-hole symmetry, present at ev- 
ery point in the Brillouin zone. The fact that it is the 
particle-hole symmetry rather than the lack of scattering 
that protects the Dirac node is clearly revealed in the re- 
lated problem of a Schrodinger electron in the presence of 
a single Aharonov-Bohm half-flux, where the density of 
states acquires a 5 function depletion at k = Q , thus 
shifting part of the spectral weight to infinity due to re- 
markably strong scattering. The authors of Rcf. |^ also 
corrected Ref. [|2| by showing that the "lines of nodes" 
must actually be the "lines of near nodes" since true ze- 
roes of the energy away from Dirac node are prohibited 
on symmetry grounds. Still, these "lines" will act as true 
nodes in all realistic circumstances, due to extraordinar- 
ily small excitation energies. 

Marinelli et al. also showed that, if the particle- hole 
symmetry is broken, for example by introducing a non- 
Bravais vortex lattice with broken inversion symmetry 
and four vortices in the unit cell, then true lines of nodes 
can develop for values of anisotropy ratio starting already 
at ~ 5. They concluded, that the density of states is 
finite at zero energy and the semiclassical results of Kop- 
nin and Volovik [e3 might apply down to zero energy. 



For a non-Bravais lattice with two vortices per unit cell 
they found that the quasiparticle spectrum can become 
gapped. 

Very recently Ye ||2^ discussed transport properties 
of the quasiparticles described by the Dirac Hamilto- 
nian of Rcf. [|l2| and pointed out some intriguing effects 
that may take place in random vortex arrays. Also, Alt- 
land, Simons and Zirnbauer investigated general proper- 
ties of disordered Dirac operators, including vortex dis- 
order [p6|| . 

In this paper we extend the original analysis which 
was based solely on the continuum description by intro- 
ducing a tight binding "regularization" of the full lat- 
tice BdG Hamiltonian, to which we then apply the FT 
gauge transformation. Our motivation is twofold: First, 
we have found by explicit numerical computations that 
different choices of singular gauge transformation result 
in spectra which, while rather similar, are not the same. 
Within our numerical accuracy we could not tell whether 
the spectra have a very slow convergence to the same final 
result or whether they actually converge to a different an- 
swer. This will be discussed again shortly. This problem 
appears to be a conspiracy between the strong Aharonov- 
Bohm scattering from magnetic half-fluxes which tends 
to push some states of the unperturbed Hamiltonian to 
very high energies and the unbounded nature of the Dirac 
spectrum. It is an interesting issue for future study how 
to devise the cutoff in the reciprocal lattice sums of the 
linearized problem which is tailor-made for a particular 
choice of the singular gauge transformation. In this pa- 
per, we circumvent this problem entirely by regularizing 
the original Hamiltonian on a square lattice. The tight 
binding formulation regularizes the strong mixing of the 
basis vectors through the introduction of an upper and a 
lower bound on the spectrum, thus prohibiting the shift 
of the spectral weight to infinity |^^. This immediately 
solves our problem: different choices of singular gauge 
transformation now rapidly converge to identical spec- 
tra, as they should. The low energy part of the spectrum 
compares best with the original FT transformation [T^ ] 
of the linearized Hamiltonian, which might have been ex- 
pected based on its having the smoothest relative phase 
between particles and holes. 

Second, the lattice formulation allows us to study 
what, if any, role is played by internodal scattering which 
is simply not a part of the linearized description. We find 
that under special circumstances, when there is a high 
degree of commensurability between the ionic and vortex 
lattices, the interference between the nodes can lead to 
scattering which is surprisingly strong (~ ^/B) and might 
be observable in HTS. Such scattering is responsible for 
opening a gap at the Fermi surface even in the case a 
Bravais vortex lattice. In a typical situation, however, 
when the two lattices have a low degree of commensura- 
bility or are of different symmetry and particularly when 
weak thermal or quenched disorder is included, the in- 
ternodal scattering effectively disappears. We diagonal- 
ize the tight binding Hamiltonian numerically for various 
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order parameter symmetries and both square and trian- 
gular vortex lattices. Our treatment provides an access 
to the entire quasiparticle energy spectrum together with 
displaying the utility of the FT transformation in ana- 
lyzing gapped superconductors (e.g. s- or Px+iy- wave), 
which are a priori inaccessible through the linearization. 
We are therefore able to present a unified treatment of a 
general, both conventional and unconventional, strongly 
type-II superconducting pairing in the vortex state. 



II. BDG HAMILTONIAN AND THE SINGULAR 
GAUGE TRANSFORMATION: LOW ENERGY 
PHYSICS OF QUASIPARTICLES AND 
VORTICES 

Because of the nonlocality inherent in the supercon- 
ductors with higher angular momentum pairing, their 
Hamiltonians are most naturally formulated on a discrete 
real space lattice representing the underlying crystalline 
lattice of the compound in question. Quite generically, 
the simplest lattice Hamiltonian which allows pairing to 
occur in s-, p- and d-wave channels is the tight binding 
model with the on-site or nearest neighbor attraction be- 
tween electrons. Conventional mean field Hartree-Fock- 
Bogoliubov decoupling of the interaction term then leads 
to the BCS-type lattice Hamiltonian of the form 



H = 



h 

A* 



A 

-h* 



(2) 



where 



and 
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s 



(4) 



The sums are over nearest neighbors and on the square 
lattice S = ±x,±y; A(r) is the vector potential associ- 
ated with the external magnetic field B, e^? is Fermi en- 
ergy, and ss is an operator which is defined by its action 
on a general function u{r) so that ss u{r) — u{r + d). The 
operator fjs depends on the type of pairing as discussed 
later. 

Quasiparticle wavefunction is a rank two spinor in the 
Nambu space, i^'^ir) = [^(r), t;(r)], and obeys the BdG 
equation 



Hip — ^ip- 



(5) 



Besides relying on conventional mean-field BCS decou- 
pling, Hamiltonian contains two additional approxi- 
mations. First, we have assumed that the order parame- 
ter magnitude is constant and equal to Aq everywhere in 
space. This is essentially the London limit llj] which is 



expected to be valid in the regime of low fields, B <C i?c2, 
when vortex cores comprise negligible fraction of the sam- 
ple. Second, we approximated the phase of the order pa- 
rameter </)a(r), which is a nonlocal field associated with 
a bond between two neighbor sites, by the average of the 
phases associated with the attached sites. 



</>5(r)^i[0(r) + </)(r + <5)]. 



(6) 



This replacement is discussed in more detail in the Ap- 
pendix A and we expect it to be very accurate far away 
from the vortex cores where the phase varies slowly, but 
inadequate in the core. Hamiltonian is therefore use- 
ful when considering quasiparticle properties in a dilute 
vortex lattice, which is the main focus of this work. To 
study properties of the core region one must explicitly 
treat the order parameter amplitude variation and non- 
locality of its phase as done e.g. in Refs. |l^,|l^. Surpris- 
ingly, however, we shall see below that even the present 
approximation yields results for the core region that are 
qualitatively correct. 



A. Continuum formulation 

In many cases our main interest is directed at the long 
wavelength and low energy or low temperature proper- 
ties. It is precisely in this respect that the quasiparticle 
excitations in an unconventional superconductor differ 
most dramatically from its s-wave counterpart. Under 
these circumstances it is desirable to consider a contin- 
uum version of the BdG Hamiltonian. For a d-wave su- 
perconductor such a continuum Hamiltonian was derived 
by Simon and Lee |27| . It turns out, however, that as 
written in Ref. [ p7| this Hamiltonian is not gauge invari- 
ant |2^. At fault is the off-diagonal term representing 
the c?-wave pairing operator, which does not transform 
properly under the U(l) gauge group. In Appendix A 
we have derived the gauge invariant form of this pairing 
operator for a pure (1^2 _y2 superconductor and have out- 
lined how such a derivation can be carried out for other 
unconventional pairing states. The continuum Hamilto- 
nian reads: 



n = 



He 

A* 



A 

-nr. 



(7) 



with Tie = 27^(P^ f A)^ — ei? and p = —ihV the momen- 
tum operator. If we follow the convention and choose the 
coordinate axes in the direction of gap nodes the gauge 
invariant d-wave pairing operator has the form 



A A(r)}} -I- -p/A(r)(|5,p,», (8) 

where p_f is the Fermi momentum, is the phase of the 
superconducting gap A(r) and curly brackets represent 
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symmetrization, {a,b} — ^{ab + ba). The above pair- 
ing operator resembles the famihar Simon-Lee form ex- 
cept for the last term which is necessary to restore the 
full gauge invariance. We emphasize that expression (||) 
is valid for uniform gap amplitude; otherwise additional 
terms which involve derivatives of the amplitude appear. 

We now use this Hamiltonian as the starting point in 
our discussion of low energy quasiparticles in the presence 
of magnetic field. Operationally, the main difhculty en- 
countered when solving for eigenstates of (0) in the vortex 
state is the nontrivial structure of the order parameter 
phase field (t>(r), which is constrained by topology to wind 
by 2tt around the center of each vortex. Ideally, we would 
want to get rid of this phase to make the problem look 
as close as possible to the reference solution in which the 
phase can simply be set to zero. If (t>(r) is a pure gauge, 
i.e. V X V0(r) = 0, this is easily accomplished by per- 
forming a gauge transformation 



gi0(r)/2 







,-i0(r)/2 



(9) 



After this transformation the BdG Hamiltonian becomes 



PxPy 



where 



, . 1 , e . , 
Vs(r)--(-V(/.--A) 
m 2 c 



(10) 



is the conventional superfluid velocity. We recognize the 
term containing V0 as the Doppler shift of quasiparticles 
in a locally uniform superflow |^,||. 

However, if (j)(r) contains vortices the situation is far 
more subtle: the vector field V(/)(r), while still locally 
uniform, acquires a global curvature, i.e. 



V X V<^(r) = 27rz ^ S{r - r,) ^ 0, 



(11) 



where {vi} denotes vortex positions. Consequently, in 
the vortex state it is no longer possible to eliminate the 
superconducting phase by the above transformation (|^) 
and obtain a Hamiltonian describing quasiparticles cou- 
pled to the locally uniform superflow. Formally this can 
be seen from the fact that in the presence of vortices 
transformation is not single valued. In principle such 
multiple valuedness of the resulting Hamiltonian could 
be handled by introducing compensating branch cuts in 
the quasiparticle wavefunctions. In practice, however, it 
is far more desirable to avoid any such complications in 
the flrst place. 

We follow FT and perform a "bipartite" singular 
gauge transformation. 









'i4>h(r) 



(12) 




magnetic unit cell 



,-o. 
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(a) (b) 

FIG. 1. Example of A and B sublattices for the square (a) 
and triangular (b) vortex lattice. 



where <^e(r) and (fihi^) are two functions satisfying 

<^e(r)+0„(r) = 0(r). (13) 

This more general transformation also eliminates the 
phase of the order parameter from the pairing term of 
the Hamiltonian but 0e (r) and 0^ (r) now can be chosen 
in a way that avoids multiple valuedness and the asso- 
ciated complications. The way to accomplish this is to 
assign the singular part of the phase field generated by 
any given vortex to either (t>e(j) or (phij), but not both 
as is done by symmetric transformation (^). Physically, 
a vortex assigned to </'e(r) will be seen by electrons and 
be invisible to holes, while vortex assigned to (j)h{'r) will 
be seen holes and be invisible to electrons. 

In practice we implement the above transformation by 
dividing vortices into two groups A and B, positioned at 
{rf} and {rf } respectively (see Fig. We then define 
two phase fields 0a (r) and 0_B(r) such that 



V X V0^(r) 



2^z^5(r-rf), 



^i = A,B, (14) 



and identify 0e = (j)A and 4>h = (j^B- Comparison with Eq. 
( |Tl| ) confirms that this choice of 0e(r) and 0/i(r) satis- 
fies the condition (p^), possibly up to some unimportant 
nonsingular phase which can be always transformed away 
by the conventional gauge transformation (|9[) . 
The transformed Hamiltonian becomes lEyl 



mv 
I) 



A\2 



D 



HP- 



,B\2 



2m 



with^ = A^[p, + f(vf,-vB)][p^ + f{vA_yB^)^ + ^^ , 



y) and 



-(W(/)p 
m 



A.B. 



(15) 



From the perspective of quasiparticles and vf can 
be thought of as effective vector potentials acting on 
electrons and holes respectively. Corresponding effec- 
tive magnetic field seen by the quasiparticles is B^^^j = 
— — (V X v^), and can be expressed using Eqs. ( p^ ) and 
(gas 
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B 



A,B, 



(16) 



sy 



"syJ 



(18) 



where B = V x A is the physical magnetic field and 
(/)o — hc/e is the flux quantum. We observe that quasi- 
electrons and quasi-holes propagate in the effective field 
which consists of (almost) uniform physical magnetic 
field B and an array of opposing delta function "spikes" 
of unit fluxes associated with vortex singularities. The 
latter are different for electrons and holes. As discussed 
in ||l2| it is desirable to choose A and B vortices in such a 
way that the effective magnetic field vanishes on average, 
i.e. (Bgjj) = 0. This translates to a simple requirement 
that we have precisely one flux spike (of A and B type) 
per flux quantum of the physical magnetic field. In that 
case flux quantization guarantees that the right hand side 
of Eq. ( |l6|) vanishes when averaged over a vortex lattice 
unit cell containing two physical vortices. It also implies 
that there must be equal numbers of A and B vortices in 
the system. 

The essential advantage of the choice with vanishing 
(B^jj) is that and vf can be chosen periodic in space 
with periodicity of the magnetic unit cell containing one 
electronic flux quantum hc/e. Notice that vector poten- 
tial of a field that does not vanish on average can never 
be periodic in space. Condition (B^fg) = is therefore 
crucial in this respect. 

The singular gauge transformation (|2|) maps the origi- 
nal Hamiltonian of fermionic quasiparticles in finite mag- 
netic field onto a new Hamiltonian which is formally in 
zero average field and has no singular phase winding in 
the off-diagonal components. The main advantage of the 
FT transformation is that it eliminates the need to in- 
troduce branch cuts into fermionic wavefunctions: these 
wavefunctions remain single valued while the physical ef- 
fect of the branch cuts is now promoted to the level of 
the Hamiltonian, where it is represented by a statistical 
gauge potential describing a half-flux Aharonov-Bohm 
scattering. This situation bears some similarity to the 
fractional quantum Hall effect (FQHE). Here, the com- 
posite fermion [^,0 is created by attaching a flux tube 
to the electron. The details, however, are quite different. 
In the present case it is the superconducting condensate 
that creates the fictitious "flux spikes" which then on 
average exactly neutralize the physical applied magnetic 
field. Unlike in FQHE, the fiuxes are stationary and we 
are generally in the limit where there is a large number 
of electrons per flux. 

To facilitate further insights into the physics of the low- 
energy quasiparticles we now linearize the transformed 
Hamiltonian in the vicinity of one of the four nodes of 
the gap function on the Fermi surface. Following Simon 
and Lee [pTl we obtain H n — Ti-o + H' , where 



^ ^ / VfPx VAPy 

is the free Dirac Hamiltonian and 



(17) 



Here vp is the Fermi velocity and wa = Aq/pf denotes 
the slope of the gap at the node. 

T-Ln can be viewed as a relativistic Hamiltonian for 
a 2-1-1 massless "Dirac" fermion and can be rewritten 
accordingly as 

■Hn = vf{Px + ax)Tz + VA{Py + ay)Ti + mvFVsx , (19) 

where are Pauli matrices, — ^(v^ -|- v^) = 
^(|V(?!)— I A) is the conventional superfluid velocity and 
a = ^(vf ^ vf ) = |(V(/)A — V(/)_b) is the internal gauge 
field. We observe that Vs couples to the Dirac fermions 
as a scalar potential while a couples as a vector poten- 
tial. The Dirac "magnetic field" b = V x a produced 
by this vector potential is highly unusual: it consists of 
delta function spikes located at the vortex centers and it 
vanishes on average when the numbers of A and B vor- 
tices are equal. Each spike carries precisely one half of 
the conventional electronic flux quantum 0o there- 
fore, although comprising a set of measure zero in the 
real space, the flux spikes lead to maximal Aharonov- 
Bohm scattering and have strong effect on the quasipar- 
ticle spectra. In particular, note that the cyclotron mo- 
tion in a Dirac cone arises entirely through b = V x a 
and does not include explicitly the actual magnetic field 
B = V X A. Such half-flux scattering is time-reversal 
invariant and cannot lead to Dirac- like (or any!) Landau 
level quantization. 



B. Internal gauge symmetry 

Spectral properties of the continuum linearized Hamil- 
tonian ( p"7| - |l^ ) have been analyzed in great detail [Il2|j22[ 
and initial investigation of its transport properties has 
been presented |25|. Here we wish to point out a pe- 
culiar property of the linearized Hamiltonian as regards 
the choice of A and B subsets of vortices, which seems 
to have been overlooked thus far. 

Logic dictates that all measurable quantities must be 
independent of our choice of A and B. This is because 
there should be no physical distinction between A and B 
vortices, the assignment being completely arbitrary. The 
freedom of assignment of vortices into A and B subsets 
represents an internal gauge symmetry of the problem 
closely related to the fact that electrons condense in pairs 
and therefore vortices carry half of the electronic flux 
quantum hc/e. 

To explicitly test this internal gauge symmetry we have 
diagonalized the linearized Hamiltonian ( [l7[{18| ) using the 
Bloch wave technique described in Ref. |12[ for two dis- 
tinct choices of A-B sublattices as illustrated in Fig. |^. 
We used a unit cell containing 4 vortices in order to be 
able to compare the band structures for the two cases 
directly on the same Brillouin zone. The corresponding 
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FIG. 2. Two sublattices, "ABAB" and "AABB" used to 
investigate the internal gauge symmetry of the FT transfor- 
mation. The shaded region marks the unit cell used in the 
numerical diagonalization. 



band structures for the isotropic case (ao = vf/va = 1) 
are displayed in Fig. ||. We observe that although quali- 
tatively similar, their detailed features are different and 
so are the associated densities of states. Similar situa- 
tion occurs for other values of Dirac cone anisotropy a d 
and other symmetries of the vortex lattice, although the 
case shown in Fig. ^ is an extreme example of the differ- 
ences. This is a surprising and unexpected result whose 
ramifications we do not fully understand at the present 
time. 

We expended considerable effort to verify that the dif- 
ference between the two band structures is not a triv- 
ial artifact of our diagonalization procedures. Rather, it 
appears to be associated with the pathological r^^l"^ 
behavior of the Dirac wavefunctions in the vicinity of a 
vortex center, which is presumably difficult to mimic us- 
ing finite number of Bloch waves that are used as basis 
states in our numerical diagonalization. We also note 
that the Aharonov-Bohm scattering induced by the half- 
flux spikes is extraordinarily strong. As shown by Moroz 
psf , in the case of ordinary Schrodinger electron it causes 
a transfer of spectral weight from zero energy up to infi- 
nite energy. 

The problem is clearly inherent only to the linearized 
BdG Hamiltonian. In the following Section we show that 
no such problem arises in the lattice version of the BdG 
Hamiltonian. This is presumably because the spectrum 
is bounded (by the tight binding bandwidth) in this case 
and therefore all states are accounted for in the numerical 
diagonalization. Also, the lattice spacing acts as a nat- 
ural short distance cutoff which regularizes the behavior 
of the wavefunctions at the core. 

We have also solved the linearized problem by directly 
discretizing the Hamiltonian jl^p^ on a square grid in 
the real space, a technique similar to that described in 
Ref. 1^^. The problem persists in this case. We conclude 
that the problem appears to be a caused by a conspiracy 
between the strong Aharonov-Bohm scattering and the 
unbounded nature of the Dirac spectrum of Hamiltonian 
( p7[]r^ ). While we believe that there exists a regulariza- 
tion scheme which would resolve the problem within the 
linearized formulation, our attempts to construct such a 
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FIG. 3. The band structure (top) and DOS (bottom) of 
the linearized Hamiltonian with two choices of sublattices: 
ABAB (thick line) and A ABB (thin line). 

scheme were unsuccessful so far. We leave it as an inter- 
esting problem for further investigation. 



C. Lattice formulation 

It is straightforward to apply the FT singular gauge 
transformation ( |l^ to the lattice BdG Hamiltonian (||). 
One obtains Hamiltonian Hn of the form 



where 



V,^(r) - ~ ^A) • dl , ix^A,B, 



and 



(20) 



(21) 



(22) 



We notice that the integrand of Eq. ( pi[) is proportional 
to the superfluid velocities defined by Eq. (|l5|) in 
connection with the continuum Hamiltonian. Appendix 
B describes an efficient method for calculation of these 
quantities in the vortex lattice. 

The main benefit of re-framing the original problem in 
this way is the explicit gauge invariance. In the case of a 
periodic arrangement of vortices the Hamiltonian is pe- 
riodic with periodicity of a magnetic unit cell containing 
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a pair of A and B vortices. In what follows we consider 
square and triangular vortex lattices with two physical 
vortices per unit cell as illustrated in Fig. |l|. The vortex 
center is always placed at the center of the plaquette of 
the underlying tight binding lattice. 

Following FT we use the familiar Bloch states as the 
natural basis for finding the eigenvalues of TIm specified 
above. In particular we seek the eigensolution of the BdG 
equation Ti-Nip — £^ in the Bloch form 



$«k(r) 



C/nk(r) 
Kk(r) 



(23) 



where {Unk, Kik) are periodic on the corresponding unit 
cell, n is a band index and k is a wave vector from the 
first Brillouin zone. Bloch wavcfunction <I'„k(r) satisfies 
the "off-diagonal" Bloch equation 



"^k^^nk ^nk*^nki 



-ik-T' 



ik-r 



with the Hamiltonian of the form TYk = 
In the following Section we describe specific forms of such 
Hamiltonians for the s-, p- and d-wave symmetries and 
discuss their solutions. 



III. NUMERICAL RESULTS 

A. s-wave pairing 

In the case of s-wave pairing the operator 775 takes the 
form fjs = \, and the Hamiltonian simplifies consider- 
ably. In particular, the off-diagonal terms become simply 
Ao, and Hn reads 



ss + ep 



(24) 



It is interesting to note that in the limit of high quasi- 
particle energy, e ^ Ao, the off-diagonal terms become 
irrelevant, and the equations for the electron and hole 
part of the Nambu wavefunction decouple. We recover 
a Hamiltonian describing holes and electrons in a uni- 
form magnetic field pierced by a lattice of counteracting 
full Aharonov-Bohm magnetic flux tubes with unit flux 
quanta hc/e concentrated at the set of point cores. The 
solution is just the familiar Schrodinger Landau levels 
(not to be confused with Eq. |l]) because the full elec- 
tronic flux has no effect on the particle energy spectrum 
p2[ . This result is expected from the outset since at 
high energies the quasiparticles behave as normal elec- 
trons or holes, which know little about the condensate. 
These high energy quasiparticles experience effectively a 
uniform magnetic field and move along cyclotron orbits. 
Similar argument holds for any pairing symmetry and 
we expect Landau level quantization of the quasiparticle 
spectrum at energies much larger than Ao. 
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FIG. 4. Magnetic Brillouin zone for the square vortex lat- 
tice with the corresponding notations used in the discussion 
of the quasiparticle band structure. 

We have numerically diagonalized the above Hamilto- 
nian making use of the standard LAPACK diagonaliza- 
tion routine. We considered a tight-binding lattice of 
10 X 10 sites, which turns out to be sufficiently large to 
analyze the CdGM regime. The corresponding magnetic 
field B = 1/(100^^) in units of unit flux hc/e per unit 
area, the superconducting gap Ao = t and the chemical 
potential ep = —2.2t, assuring an approximately cylin- 
drical Fermi surface. The resulting spectrum for the Bril- 
louin zone displayed in Fig. ^ and density of states for the 
square vortex lattice are shown in the Fig. ^ The B = 
spectrum has the usual BCS form with a full gap Aq . The 
additional features at 2.1Ao and 2.4Ao are remnants of 
the band edge and the van Hove singularity respectively 
present in the normal state spectrum. 

The magnetic field induces low-energy states within 
the gap, which become localized in the vortex cores. 
These are CdGM states |l^ dispersed into bands. At 
low energies, the bands are very narrow signaling strong 
concentration of the wavefunctions at the vortex cores 
and insignificant overlaps among the states at neighbor- 
ing vortices. This fact justifies the chosen parameters. At 
energies less, but comparable to Aq, the bands are broad- 
ened due to increasing overlap among the wavefunctions. 
In is interesting to note that CdGM bound states ap- 
pear despite the fact that our model assumes constant 
order parameter amplitude and the effective core size is 
the tight binding lattice spacing 6. The small size of 
the core causes the lowest bound state to be pushed to 
rather high energy and also that only few bound states 
can be resolved with our numerical accuracy. For energies 
e 3> Ao the spectrum exhibits Landau level quantization, 
as expected from the argument presented above. 

It is appropriate to illustrate the pitfall lurking in guise 
of the symmetric transformation widely used in the litera- 
ture. At the first glance, perhaps the most natural choice 
for removing the phases from the off-diagonal terms is 
setting in Eq. (H) (^e(r) = (j)h{r) = (t){r)/2. Note that in 
this case the transformation 



U 



(,-i'l'{r)/2 



(25) 



is not single valued and neither are the resulting wave- 
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FIG. 5. Top: Quasiparticle band spectrum for an s-wave 
superconductor in the presence of the external magnetic field 
B=l/1005^, and Ao = t, ep ~ —2.2t. Bottom: correspond- 
ing DOS. Note the bound Caroli-Matricon bands at energies 
below the gap and the Landau levels at energies e ^ Aq. 



functions. Nevertheless, ignoring these facts, the Hamil- 
tonian Ti. pf becomes 



iVsir 

Ao 



ss - ep 



Ao 



with 



(iv0-|;A).dl 



(26) 



(27) 



In the hmit of high quasiparticle energies the equations 
again decouple, but now they describe a quasiparticle 
moving in a uniform magnetic field pierced by half elec- 
tronic flux quanta hc/2e canceling the overall field. These 
half-fiuxes cause significant Aharonov-Bohm scattering 
and cannot be ignored. As shown in Ref. |^2|, the spec- 
trum for this problem is not that of Landau levels; there 
is a significant dispersion. We see that symmetric trans- 
formation leads to the results that are manifestly incor- 
rect. Again, this argument is independent of the pairing 
symmetry. 



B. 



p-wave pairing 



We follow Matsumoto and Sigrist |g] and for simplic- 
ity assume that the prototype p-wave superconductor is 
two dimensional and has a cylindrical Fermi surface. We 
further assume for simplicity that it is strongly type-II, 
although Sr2Ru04 is not of this type. We restrict the 
Hamiltonian to one of the two degenerate states, Px+iy 
and ignore the Px-iy part. For p-wave {px + ipy) we have 



±ss 



if (5 
if (5 



±x 



(28) 



where ss u{v) — u{y + 8). The Hamiltonian becomes 



-tEs 
AoE. 



ss - ep 



AoE. 



(29) 



where the phases V^(r) are defined by Eq. 



Asir) = \jj (V0A - V<^b) • dl 



and 



(30) 



We chose ep = —2.2t which yields approximately cir- 
cular Fermi surface with the superconducting gap, to a 
good accuracy, uniform everywhere on the Fermi surface. 
As in the case of s-wave, the value of Ao is set equal to t. 
The resulting spectrum and density of states are shown 
in the Fig. (^. 

The spectrum again reveals bound vortex states broad- 
ened into a band. In contrast to the s-wave case we now 
have a state at zero energy. These results are what is 
expected on the basis of our understanding of a single 
^»-wave vortex MM . 
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FIG. 6. Top: Quasiparticle band spectrum for a pa;+iy-wave 
superconductor in the presence of the external magnetic field 
B=l/100(5^, and Ao = t, ep = ~2.2t. Bottom: The corre- 
sponding DOS. 



10 20 30 40 50 60 70 
FIG. 7. Local density of states at i5 = for the 
d^2_y2-wa,ve superconductor with square arrangement of vor- 
tices. The plot is in units of the tight binding lattice con- 
stant 5. Bright regions represent maxima while the dark 
regions represent minima. The parameters are ep = 0, 
an = t/ Aq — 4. The strong overlaps of the low energy quasi- 
particle wave functions along the four nodal directions cause 
appreciable interference effects which in turn influence the 
character of the quasiparticle spectrum. 



C. d-wave pairing 

To model the high-temperature superconductors such 
as YBa2Cu307, we assume that coupling between the 
Cu-0 planes is weak and to the leading approximation 
can be ignored. On the tight binding lattice the d-wave 
pairing is given by A = 2Ao[cos(/ca;(S) — cos{ky5)] which 
determines the form of the lattice operator to be: 



Vs 



sg if 5 = ±x 
-ss if (5 = ±y 



(31) 



where as before ss u{r) ^ u{r + 5). With this definition 
of ris the Hamiltonian for d-wave pairing has the same 
form as Eq. (p9|). 

The results presented in this section correspond to 
the magnetic field 0o/(16OO(5^) for square vortex lat- 
tice and (/)o/ (1500(5^) for triangular vortex lattice where 
00 = hc/e = 4.137 x 10^ tA^ and S is the tight binding 
lattice constant. Taking (5 = 4 A, as in YBa2Cu307, this 
corresponds to physical field of 16 T. The above param- 
eters were chosen for computational efficiency, but we 
did not see any qualitative difference down to the fields 
as low as i^o/ (4900(5^) corresponding to a magnetic unit 
cell of 70S X 706 and a field of 5.2 T. Numerical diag- 
onalization was performed using the ARPACK package 
routines for sparse matrices. This algorithm provides a 
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FIG. 8. Low energy part of the quasiparticle band spec- 
trum for the square vortex lattice. The parameters are ef = 0, 
an = 4. Top: example of a gapless spectrum for I = 385. 
Note that the node at Q is moved away from the original 
Dirac F point. This effect is a result of a uniform gauge 
"boost" associated with the choice of vortex unit cell and dis- 
appears for a unit cell with four vortices. Bottom: example 
of a gapped spectrum with I = AOS. 



set of low-lying eigenvalues and allows handling much 
larger systems than the full diagonalization used in s- 
and p-wave cases. 

We find that the quasiparticle wave functions exhibit 
significant dependence on the symmetry of the vortex 
lattice. For the square lattice, the overlap among the 
wave functions corresponding to different nodes is ap- 
preciable and there are strong interference effects along 
the I a; I = \y\ diagonals, i.e. the directions in the real 
space where A(k) vanishes. This is illustrated in Fig. |^. 
For certain commensurability of the tight binding and 
the square vortex lattices, the interference effects are re- 
sponsible for opening a gap at Fermi energy, while for 
the complementary set of lattices at different commen- 
surability factors, the spectra are gapless at the Dirac F 
point. 
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FIG. 9. Low energy part of the quasiparticle density 
of states for an d^2_y2-wave superconductor with square ar- 
rangement of vortices for two different interference cases in 
which \^F ■ d = 2nn (solid line) and kF • d = (2n -I- l)7r 
(dashed line), n being an integer, kF = (-|,|-)a Fermi vector 
at nodal points, and d is the primitive vortex lattice vector. 
Notice the appearance of the gap in the density of states for 
kF ■ d = 2n7r. Plotted on arbitrary scale, energy is in units of 
t. The parameters are ep = 0, an = 4. 



Figs. H and|^ show the low-energy band structures and 
the low-energy density of states for the square lattice. 
The two system sizes shown illustrate the commensura- 
bility effect: if the scalar product between the Fermi vec- 
tor along the nodal direction ki? and the vortex primitive 
Bravais lattice vector d is an even integer times tt, the 
spectrum develops a gap, while it remains gapless if this 
product is an odd integer times tt. The same effect is seen 
at higher Dirac anisotropy ao = t/Ao = 10 (Figs. |l^, |l^) 
and ai) = 15 (Figs. ^ 

These interference effects persist down to low magnetic 
fields where the interference gaps scale as ^ \/B (see 
the next section and Fig. |l^). In the case of a triangu- 
lar lattice, the interference effects were greatly reduced 
(Fig. |l4[ ) and no commensurability dependence was ob- 
served. We find the spectrum to be gapless at half filling 
in this case (see Fig. [l5| ). 

Finally we note that we have explicitly verified that 
identical spectra (to within numerical accuracy) are 
found irrespective of our assignment of the A-B sublat- 
tices. This finding confirms that the choice of A-B vor- 
tices is an internal gauge symmetry of the problem, as 
one would expect on general grounds. 



IV. DISCUSSION 
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FIG. 10. Low energy part of the quasiparticle band spec- 
trum for the square vortex lattice. The parameters are ef = 0, 
QfD ~ 10. Top: gapless spectrum for I = 385. Note the in- 
crease of the dispersion in the QM direction with increase 
of the Dirac anisotropy q_d. Bottom; gapped spectrum for 
I = 405. 

A. Comparison of continuum and lattice results 

The results of Sections II and III show that, un- 
der generic conditions, the spectrum of linearized Dirac 
Hamiltonian provides a reasonable approximation to the 
low energy part of the spectrum of the full BdG Hamilto- 
nian. The Dirac nodes are preserved provided that there 
is no commensuration between kj^ at the node and the 
primitive vector of the vortex lattice d, and thus the in- 
ternodal scattering can be neglected. The overall shape 
of the energy bands is also qualitatively and quantita- 
tively similar for the two cases. 

Previous investigations of the linearized Hamiltonian 
|^ , p2| established that the spectrum becomes quasi one 
dimensional and lines of nodes appear |33[ in the Brillouin 
zone for large Dirac cone anisotropy ao > 14, leading to 
finite DOS at the Fermi level. Inspection of Fig. |l^ sug- 
gests that similar effect takes place in the full BdG Hamil- 



Energy [t] 

FIG. 11. Low energy part of the quasiparticle den- 
sity of states for an ci^2_y2-wave superconductor with square 
arrangement of vortices for two different interference cases: 
I = 38(5 (dashed) and I = 405 (solid). Plotted on arbitrary 
scale, energy is in units of t. The parameters are ep = 0, 
ao = 10. 

Ionian, although the one dimensionality is somewhat less 
pronounced and is restricted to the immediate vicinity 
of the node. Furthermore, the lines of nodes never quite 
form (although the tendency is clearly visible along the 
line r — > M) and the DOS remains zero at the Fermi 
level. In the above discussion one needs to bear in mind 
that the vortex lattice considered here has been rotated 
by 45° relative to Ref. H. 

B. Scaling of the Energy Spectrum 

The wavefunction interference effects among the four 
nodes, which are responsible for opening of the "inter- 
ference" gaps visible in some of our spectra, seem to be 
in contrast with the results obtained for the linearized 
Hamiltonian by FT |12| and more recently by Marinelli, 
Halperin, and Simon |23|, where any internodal interac- 
tion is ignored and assumed insignificant. Furthermore, 
Marinelli et al. advanced strong analytic arguments that 
in the presence of particle-hole symmetry the linearized 
Hamiltonian retains the Dirac node at the F point and 
does not develop a gap to order O (^~^), where I is the 
magnetic length. We found that the "interference" gaps 
in the quasiparticle spectrum in the case of the square 
lattice scale with magnetic field as ^fB ~ This is 

shown in Fig. |l6|. The reason for this can be understood 
from the scaling of the wavefunctions of the linearized 
Hamiltonian. We find that there is a r~^/^ divergence in 
the asymptotic solution of the wavefunction around one 
vortex. This strong concentration of the wavefunction 
around the vortex makes the contribution from the term 
quadratic in the supcrfluid velocity particularly enhanced 
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FIG. 12. Low energy part of the quasiparticle band spec- 
trum for the square vortex lattice. The parameters are ef = 0, 
an ~ 15. Top: gapless spectrum for I = 385. Note the in- 
crease of the dispersion in the QM direction with increase 
of the Dirac anisotropy q_d. Bottom; gapped spectrum for 
I = 405. 

and, independent of regularizing the wavefunctions to 
eliminate the divergences at the core, the contribution 
to the gap is significant. We can then extract the depen- 
dence of the wavefunction on the magnetic length I just 
from dimensional analysis. The wavefunction must have 
units of length"^ therefore ~ (rZ)~^/^. One can see, 
that the matrix element, and consequently the gaps, will 
in general scale as for the terms beyond the linearized 
Hamiltonian. This dependence is extremely difficult to 
obtain from the plane wave expansion of the wavefunc- 
tions. 

Our numerical results strongly suggest that there is a 
characteristic oscillation in the gap of the spectrum de- 
pending on the commensurability of the magnetic lattice 
and the underlying ionic lattice. This can be interpreted 
as the inter-nodal scattering. The interaction between 
the quasiparticles at different nodes is responsible for 
opening the gaps at the Fermi surface. The effect of the 
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FIG. 13. Low energy part of the quasiparticle den- 
sity of states for an ci^2_y2-wave superconductor with square 
arrangement of vortices for two different interference cases: 
I = 385 (dashed) and / = 405 (solid). Plotted on arbitrary 
scale, energy is in units of t. The parameters are ep = 0, 
an — 15. 
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FIG. 14. Displayed is a typical low energy modulus of a 
quasiparticle wavefunction for the d^2_y2-wa.ve superconduc- 
tor with triangular arrangement of vortices. The plot is in 
units of the tight binding lattice constant 5. Bright regions 
represent maxima while the dark regions represent minima. 
The parameters are ep = 0, an = 4. This plot illustrates the 
reduction of the interference effects for triangular vortex lat- 
tice. As a consequence the gap in the quasiparticle spectrum 
does not emerge as shown in Fig. 
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FIG. 15. Low energy part of the quasiparticle density of 
states for an d-wave superconductor with triangular arrange- 
ment of vortices. Plotted on arbitrary scale, energy is in units 
of £. The parameters are ep = 0, ao = 4. 



intra-nodal scattering on these gaps on the Fermi surface 
is, however, absent since for certain commensurability of 
the magnetic and ionic lattices there is no gap. Thus we 
conclude that the effect is purely due to the internodal 
scattering mediated by the terms beyond the linearized 
Hamiltonian. The sensitivity of the gaps to the commen- 
surability of the ionic and magnetic lattices is supported 
by the results with triangular vortex lattice, in which case 
the spectrum remains gapless as there is no commensu- 
rability between ionic and vortex lattice. This supports 
the view that the internodal scatting alone is responsible 
for the presence of the gap at the Fermi surface. 
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FIG. 16. The magnitude of the interference gaps vs. mag- 
netic length I exhibits scaling. Ao = 0.25i, ep ~ 0. 



the Doppler shift terms and reads 



SS + EF 



(32) 



The density of states obtained by diagonalizing Tis 
is shown Fig. 
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It is significantly different from the 
results presented in Section III. This clearly demonstrates 
that there is an essential piece of physics missing from 
the Hamiltonian which contains only the Doppler shift 
effect, and consequently from the frequently encountered 
semiclassical approximation to such Hamiltonian W. 



V. CONCLUSIONS 



C. Comparison with the Doppler-Shift-Only Results 
for d-wave gap 

One of the key insights gained from the FT trans- 
formation is that the familiar and often used Doppler 
shift approximation is not sufficient to describe the 
quasiparticle dynamics in the vortex lattice. While the 
Doppler shift enters at the level of a linearized Dirac 
Hamiltonian as a periodic scalar potential, there is also 
an effective wecior potential a (Sec. II), which originates 
from the global curvature of the superflow in the pres- 
ence of vortices. This new vector potential term leads 
to additional strong magnetic half-flux scattering across 
the whole energy spectrum. It is instructive to compare 
our results with those obtained by performing the sym- 
metric gauge transformation specified by Eq. (|2^). As 
already pointed out the symmetric transformation is not 
single- valued and the resulting transformed Hamiltonian 
must be accompanied by the branch cuts imposed on its 
eigenf unctions. If one simply ignores these branch cuts 
altogether, the resulting Hamiltonian Ti.s contains only 



In conclusion, general utility of the singular gauge 
transformation for the calculation of the quasiparticle 
spectra in the vortex state of a general pairing symme- 
try was shown. Once the tight binding regularization is 
introduced, the spectrum can be computed in principle 
exactly using the Bloch states as the natural basis, al- 
though one is bound to resort to numerical calculation 
regardless of respecting the self-consistency. In the case 
of s- and p-wave symmetry we showed that the method 
applied to an array of vortices leads to results consistent 
with single vortex solutions. 

For d-wave pairing spectrum is also consistent with the 
single vortex solution from the point of view that the all 
the quasiparticle wavefunctions are delocalized and no 
bound states are observed. Additional insight is gained 
from the exact solution with respect to the continuum 
linearized version of the theory. For specific commensu- 
rability of the tight binding and square vortex lattice the 
internodal scattering mediated by the terms neglected in 
the linearized theory is found to be significant and of 
the same order of magnitude as the terms present in the 
linearized Hamiltonian. This is believed to be brought 
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FIG. 17. Comparison of the low energy part of the quasi- 
particle density of states for an d-wave superconductor with 
square arrangement of vortices with the DOS obtained from 
the Doppler-shift-only approximation. Plotted on arbitrary 
scale, the energy is in units of t. The parameters are = 0, 
I = 385, ao = 4. 



In particular the effect of static and dynamic disorder 
in the vortex positions on the quasiparticle spectra must 
be understood in order to make connection with the ex- 
perimental data. Another set of unresolved issues arises 
in connection with the zero field superconducting state 
phase-disordered by fluctuating vortex-antivortex pairs 
|p[-pl|. One would expect that the "Berry phase" term 
arising from fiuctuating vortices would infiuence in a pro- 
found way the critical behavior of the HTS system on the 
verge of becoming a Mott insulator. 
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APPENDIX A: LATTICE AND CONTINUUM 
BDG HAMILTONIAN IN A GAUGE INVARIANT 
FORMULATION 



forth by the diverging accumulation of the Dirac wave- 
functions in the vicinity of the vortex core, consequently 
giving rise to increased significance of the terms beyond 
linearization. However, since rather special conditions 
must be met for this effect to be significant, it is sug- 
gested that introduction of any perturbing agent such as 
disorder in the position of the vortices or vortex vibra- 
tions will lead to decoherence of the matrix elements for 
the internodal scattering and subsequently to the sup- 
pression of the interference effect. It is also possible, and 
in our view likely, that in a fully selfconsistent solution 
the vortex array would spontaneously undergo a slight 
spatial deformation into an incommensurate state so as 
to avoid opening gaps in the quasiparticle spectrum. In 
this respect, at chemical potential tp — 0, the results of 
the theory regularized on the tight binding lattice agree 
with the continuum linearized version. 

We have uncovered a peculiar property of the linearized 
Dirac Hamiltonian: it appears to violate the internal 
gauge symmetry associated with the assignment of A and 
B vortices. Although the final resolution of this apparent 
contradiction awaits further research, we attribute it ten- 
tatively to the unusually strong scattering of Aharonov- 
Bohm half fluxes acting in the Dirac Hamiltonian with 
unbounded excitation spectrum. This problem does not 
occur in the full BdG Hamiltonian. We conclude that the 
original "j4i3j4i?" choice of the gauge |l^ is the one most 
representative of the actual spectrum because it results 
in smoothest possible variation of phase in the vortex 
lattice. This view is also supported by the direct com- 
parison with the spectrum obtained using the full BdG 
Hamiltonian. 

Number of intriguing issues remain to be addressed. 



In this Appendix we derive the explicit form of the 
pairing operator A ( [4|) which appears in the lattice BdG 
Hamiltonian of Eq. (0). We also derive the continuum 
limit of this operator which is used to construct the BdG 
equations of Section II. Throughout our derivation we 
pay a special attention to the preservation of local gauge 
invariance. We start with the general pairing term on a 
tight-binding square lattice: 



^ A(*,j)K(*)«(j) + ^^*(j>W] 



(Al) 



<ij> 



Here A(j,j) is a complex paring potential defined on 
the nearest neighbor bonds. Such pairing term gener- 
ically arises when t-J and related Hamiltonians, which 
are thought to represent good microscopic models of var- 
ious unconventional superconductors, are treated within 
a BCS-type pairing approxi mati on. A conventional s- 
wave case follows from Eq. ( |AlD if we replace the sum 
over nearest neighbor bonds < i,j > with the sum over 
sites i. By construction, the pairing term (Al) is invari- 
ant under gauge transformations: 



u{i) -> uii)e'''^'\ 
v{i) -> w(i)e-*xW, 
A(z,j)->A(i,j>^'^«+^'^(^), 



(A2) 



where x(«) is an arbitrary non-singular function. 

The actual form of the complex function A(i,j) = 
D{i, j) eTq){i(p{i, j)) is obtained as a self-consistent so- 
lution of the gap equation in some specific gauge. We 
denote its "amplitude" and phase by D{i,j) and ip{i,j), 
respectively. For convenience, the "amplitude" D{i,j) is 
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defined so that it aheady contains the information about 
the relative orbital state of a superconductor. For exam- 
ple, in a pure d^^-y^ superconductor and at zero field, 
D{i,j) = — (+)A for < i,j > in the x{y) direction, 
where A is a complex constant. Furthermore, for the 
purposes of this paper, we assume that the actual am- 
plitude |A(i,j)| can be well approximated by a uniform 
(real) constant Aq, independent of < i,j >. This as- 
sumption is valid in the space between vortices but it 
clearly breaks down inside a vortex core. At low fields, 
Hci <C <SC Hc2, where the intervortex separation is 
much larger than the core size we expect any effect of 
the inhomogeneous amplitude to be negligibly small. 

The essential information about vortex configurations 
and self-consistent solution at a finite magnetic field is 
now stored in the bond phase ip{i,j). Near a plaquette 
containing a vortex ip(i,j) changes rapidly from bond to 
bond. Far from vortex cores, however, we expect ip{i,j) 
to be some smoothly varying function undergoing only 
small changes between neighboring bonds. Consequently, 
in the regions far away from vortex cores we can replace 
the bond phase (p{i,j) by a suitably chosen site phase 
variables (t>(i). The natural choice for the (/)(«) 's is a sim- 
ple average: 



(A3) 



where the sum over a runs over four bonds containing 
the site i. With this choice of (?!)(«) 's we can replace: 



(A4) 



Note that, given the choice of site variables (|A3|), Eq. 
(A4) is an approximation, accurate up to second order 
lattice derivatives of 0(j)'s. We could use a more elabo- 
rate representation of Lp{i,j)^s in terms of 0(i)'s so that 
(A4) is satisfied to even higher degree of accuracy. This, 
however, is entirely unnecessary in the present context, 
since our overall accuracy is precisely at the level repre- 
sented by (A3) and (A4). The replacement (A4) simpli- 
fies the pairing term of the lattice BdG Hamiltonian and 
reproduces the form of the pairingoperator A (^ used 
in our lattice Hamiltonian of Eq. (||). 

The above replacement of bond phases </3(i,j)'s with 
site phases 0(i)'s is also a necessary first step in our 
derivation of the continuum B dG Hamiltonian. Since 
both u{i) and v(i) appearing in ( |Al| ) are site fields we ex- 
pect that the continuum pairing term in unconventional 
superconductors will involve u(r) and v(r) acted upon 
by some local operator. To determine the explicit form 
of this operator we first combine (Al) and (A3) into: 



X [?i*(i)w(i + (t) + u*{i + cr)v{i)\ + h.i 



(T)-i0(i) 



(A5) 



where we have transformed the summation over bonds 
into the summation over sites. We now use ip{i, i + a) = 



+ \4>{i + a) + 0(^^0),(|A^ where 5"^ denotes sec- 
ond order lattice derivatives which are unimportant in 
the continuum limit. Also, from now on we restrict our 
attention to the most interesting case, a pure d^2_y2 su- 
perconductor. This allows us to rewrite (A_5) as: 



A(l)[-e*"^('+*)/2-«0W/2(y*(,)„(,;_^X)^„*(j_^X)„(j)) 



_g»</.(z-x)/2-#W/2^^*^^^^^^ - X) + - x)v(l)) 

_^g^0(^-y)/2-^0(^)/2^y*(.)^,(. _ + u* {i - y)v{i))] 

+ h.c., (A6) 

where A(i) = A exp(i0(z)) and x,y are unit displace- 
ments on the square lattice. Next, we expand 

gZ0(^±x(y))/2-z0(O/2 « 1 + 1 ± x(y)) - 0(z)] 

-l[^i^±^y))-cj,{^)]' + ..., (A7) 

make a transition to continuum variables u{i) — > au{r), 
v{{) -> av{r), A{z) ^ A(r), 0(z) ^ ,/>(r), ^ 
j{d^r/a'^) and use the standard definitions of lattice 
derivatives to finally obtain the continuum version of the 
pairing term: 

^ j d^r{u*{v)A{v)[{d,+ '-{dA{v)){d,+-{dA{r))v{v)] 



-[(9. + ^(a.0(r))(9, + '-{dA{Y))u*{Y)\/\{Y)v{x 
- (x -> y)] -l-h.c. . 



(A8) 



In going from ( [Aq ) to (|A8|) one encounters some lengthy 
but straightforward algebra. We found that decomposing 
the sum over nearest neighbors in (A6) into s-, p-, and 
d-wave components relative to site i facilitates the book- 
keeping and makes the computations rather efficient. All 
the relevant derivatives up to and including second order 
are kept and accounted for. Higher order derivatives do 
not appear reflecting our origi nal starting point of the 
nearest neighbor pairing only (Al). Note that a is the 
lattice spacing in our model. 

The form of the local continuum pairing operator is 
now apparent. We can view A(r) = A exp(z(/)(r)) as rep- 
resenting the center- of-mass portion of the gap function. 
The original non- locality, arising from the relative d^^-y-^ 
character of the pairing, manifests itself through "covari- 



ant" derivatives dr 



K9r</'(r)) 



the phase of A(r). Note that (^ 



where 0(r) is precisely 
8) is explicitly invariant 
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under the continuum version of local gauge transforma- 
tions: u{r) u{r) exp{ix{r)) , v{r) w(r) exp(— ix(r)), 
A(r) A(r) exp(2jx(r)). 

The off-diagonal elements of the Hamiltonian matrix 
appearing in the continuum BdG equat ions are obtained 
by taking the functional derivatives of (A8) with respect 
to u*(r) and v{r). This results in: 

- a^d^, {d^, A(r)}} + a^dy, {dy,A{r)}} - 

-lA(r)a2[(a,2</.)-(5»], (A9) 

and its hermitian conjugate. Here we used the stan- 
dard notation: {a, &} = ^{ab + ha). In performing the 
functional derivatives we have exploited the fact that all 
spatial dependence of A(r) comes through its phase, i.e. 
i9rA(r) = iA(r)(9r0(r), in line with our previous assump- 
tions. 

While our derivation starts with a familiar model of 



the lattice d-wave superconductor (Al) and naturally de- 



scribes the dx2-y2 state in actual continuum calculations 
it is often more convenient to consider a dxy supercon- 
ductor, so that either an a; or a ?/ axis coincide with a 
particular nodal direction, as in Section II. We can obtain 
the pairing term in the continuum BdG Hamiltonia n o f 
a dxy superconductor by simply rotating our result ( AS ) 
by 45°: 



APPENDIX B: PHASE FACTORS AND 
SUPERFLUID VELOCITIES 

In this Appendix we derive expressions for superfluid 
velocities and vf which enter both continuum and 
lattice versions of the BdG Hamiltonians in consideration 
in Section II. We start by taking the curl of Eq. (p^). 



V X v^f = 



2TTh 



zV<5(r-rf)-B/</.o 



(Bl) 



where (f>o = he/ e is the flux quantum, B = V x A, and 
we have used Eq. (p^. In the intermediate field regime 
the magnetic field distribution is to an excellent approx- 
imation described by the conventional London equation 



-0oz^(5(r-rO, 



(B2) 



where A is the London penetration depth and the sum 
now runs over all vortex positions. The London equa- 
tion is easily solved by going over to the Fourier space, 
obtaining B(r) = (27r)-2/ d2fce'*' '"Bk with 



Bl 



1 + A2fc2 



(B3) 



/ d\{u*{v)A{v)[{dx + \{dx^{v)){dy + \{dy<i^{v))v{v)] 



+ [{dx + \{dx4>{v)){dy + '-{dyC^{v))u*{v)\A{v)v(v) 



+ {x y)] ^' h.i 



(AlO) 



Similarly, by taking functional derivatives we obtain 
the off-diagonal matrix elements of the continuum BdG 
Hamiltonian operator: 

- ^{dx. {dy, A(r)}} - a\dy,{dx. A(r)}} - 

-'-A{v)a'{dxdycj>) , (All) 

and its hermitian conjugate, which is precisely the ex- 
pression used in Section II, provided that we identify 
with 2a?. 

The above derivation can be easily repeated for a p- 
wave lattice Hamiltonian and is in fact only simpler. We 
therefore do not give it explicitly but trust that the d- 
wave derivation provides a sufficiently detailed prescrip- 
tion. Similarly, our derivation is straightforwardly gener- 
alized to other unconventional forms of superconducting 
pairing. 



If we now Fourier transform Eq. (Bl) we obtain 



ik X vf^ 



— ik-i 



Bk/00 



(B4) 



To solve for v^^ we take a vector product of both sides 
with zk. After substituting for Bk and some easy algebra 
we obtain 



A 2ttTi f d^k ik X z , , 



TO J (27r)2 A:2 



1 Aw + B 



2 1 + A2fc2 



^ 1 



(B5) 



and a similar expression for with Ak and _Bk inter- 
changed. Here we have defined 



A. 



E' 



-ikrf 



Eq. ( p35| ) gives an explicit formula for which can 
be evaluated for arbitrary distribution of vortices. For 
strongly type-II materials in fields well above Hd Eq. 
(B5) may be simplified further by rewriting the expres- 
sion in the brackets as 



A, 



A2fc2 



1 Bk - Aw 



1 + A2fc2 2 1 + A2fc2 



and noting that since A2fc2 ~ \^/d^ ^ 1 {d being inter- 
vortex distance), the second term can be safely neglected. 
We thus obtain 
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2'KhX'^ f (Pk ikx z 



(27r)2 1 



(B6) 



a formula used in Ref. which is valid for all practical 
purposes. Phase factors V and A entering the lattice 
Hamiltonians of Sec tion III may be obtained as simple 
line integrals of Eq. (B6). 
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